Preprint typeset in JHBP style - PAPER VERSION 



CERN-PH-TH/2010-318 



AAMQS: A non-linear QCD analysis of new HERA 
data at small-x including heavy quarks 



Javier L. Albacete^, Nestor Armesto^, Jose Guilherme Milhano'^'^, Paloma Quiroga 
Arias^ and Carlos A. Salgado^ 

^ Institut de Physique Theorique, CEA/Saclay, 91191 Gif-sur-Yvette cedex, France. URA 
2306, unite de recherche associee au CNRS. 

^ Departamento de Fisica de Particulas and IGFAE, Universidade de Santiago de 
Compostela, E- 15706 Santiago de Compostela, Spain. 

^ CENTRA, Institute Superior Tecnico (1ST), Av. Rovisco Pais, P-1049-001 Lisboa, 
Portugal. 

^ Physics Department, Theory Unit, CERN, CH-1211 Geneve 23, Switzerland . 
5 LPTHE, UPMC Univ. Paris 6 and CNRS UMR7589, Paris, France. 

E-mail addresses: javier.lopez-albacete@cea.fr, nestor.arniesto@usc.es, 
guilherme .milhanoSist .utl.pt , pquirogaSlpthe . jussieu.fr, 
carles . salgadoOusc . es. 

Abstract: We present a global analysis of available data on inclusive structure func- 
tions and reduced cross sections measured in electron-proton scattering at small values of 
Bjorken-x, x < 0.01, including the latest data from HERA on reduced cross sections. Our 
approach relies on the dipole formulation of DIS together with the use of the non-linear 
running coupling BK equation for the description of the small-x dynamics. We improve our 
previous studies by including the heavy quark (charm and beauty) contribution to the re- 
duced cross sections, and also by considering a variable flavor scheme for the running of the 
coupling. We obtain a good description of data, the fit parameters remaining stable with 
respect to our previous analyses where only light quarks were considered. The inclusion 
of the heavy quark contributions resulted in a good description of available experimental 
data for the charm component of the structure function and reduced cross section provide 
the initial transverse distribution of heavy quarks was allowed to differ (more specifically, 
to have a smaller radius) from that of the light flavors. 
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1. Introduction 

The experimental data collected in electron-proton deep inelastic scattering (DIS) exper- 
iments at small values of Bjorken-x constitute one of the most valuable sources of infor- 
mation to test and explore the high-energy limit of Quantum Chromodynamics (QCD). 
It is a well settled theoretical result that in the limit of small Bj or ken- a; or, equivalently, 
at high energies, deviations from standard collinear perturbation theory are expected on 
account of large gluon densities developing on the proton wave function. Such corrections 
can be interpreted in different ways. On the one hand, preserving unitarity of the theory 
sets an upper limit on the growth rate of the gluon densities in the proton. This limit 
is realized by the inclusion of gluon recombination processes, highly probable in a high 
density environment, into the high-energy evolution equations. Such task is best carried 
out in the framework of the Color Glass Condensate (CGC), which is equipped with a 
set of renormalization group equations, the BK-JIMWLK equations [1-8], that include the 
needed unitarity corrections. On the other hand, the interplay between radiation (linear) 
and recombination (non-linear) process gives rise to a dynamical transverse momentum 
scale, the saturation scale, which signals the onset of non-linear corrections. The presence 
of a perturbatively large transverse momentum scale would invalidate the main assumption 
of collinear perturbation theory, where partons are assumed to move collinearly with the 
hadron, with zero transverse momentum. 
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On the phenomenological side, intense activity has been carried out in order to identify 
the onset of non-hnear corrections in available data. The dipole model formulation of 
deep inelastic scattering [9, 10] is an instrumental tool in this search, since it allows for a 
relatively simple implementation of saturation effects in the description of the scattering 
process. Indeed, dipole models [11-17] successfully describe a number of features observed 
in data at small-x. Important progress has been achieved in reducing the degree of modeling 
required in phenomenological works thanks to the determination of higher order corrections 
to the BK-JIMWLK equations [18-21]. The latter provide the most complete theoretical 
tool so far to describe the small-x dynamics of the dipole scattering amplitude including 
unitarity corrections. While higher order corrections present a complicated structure and 
are not amenable to a numerical implementation, it was demonstrated in [22, 23] that 
considering only running coupling corrections to the BK equation, henceforth referred to 
as rcBK, grasps most of the higher order effects. It was also observed in these works 
that the quantitative features of the solutions of the rcBK equation, and the much milder 
x-dependence of the saturation scale [14], are compatible with the values extracted from 
purely phenomenological analysis, rising hope that the rcBK equation could be used as a 
phenomenological tool. This idea was confirmed in a previous work by some of the authors 
[24,25], where we demonstrated the ability of the rcBK equation to correctly describe the 
data available at the time on inclusive and longitudinal structure functions mcasTircd in 
e+p collisions. Later on, a similarly successful description data for the diffractive structure 
functions was presented in [26] . The rcBK equation has also been successfully employed in 
the analysis of data on total multiplicities in Au+Au and LHC Pb+Pb collisions [23, 27] 
and on single [28] and double inclusive [29] spectrum in p+p and d+Au collision performed 
at RHIC, where saturation effects arc believed to be an important dynamical ingredient. 
Thus, the rcBK equation which, we recall, follows from a strict derivation in pQCD in a 
given limit (see section 2 for details), has become the most effective phenomenological tool 
to assess the role of unitarity effects in available data in a theoretically controlled way, thus 
bridging the gap between theory and experiment. 

On the other end of the theory spectrum, DGLAP based analyses (see [30] and refer- 
ences therein) have consistently reported good fits to e+p data for ^ 1 + 4 GeV^. A 
relevant question is whether the flexibility in the initial conditions for DGLAP evolution 
is hiding some interesting QCD dynamics, namely the presence of non-linear behaviour. A 
recent work [31] from members of the NNPDF collaboration showed that there is tension 
in NLO DGLAP fits to HERA data on structure functions when data sets with below 
the estimated saturation scale of the proton were excluded from the analysis. Moreover, 
it was found that such deviations between theory and data can not be corrected by the 
inclusion of next-to-next-to leading order corrections in DGLAP evolution nor by a better 
treatment of heavy quark effects. It has been suggested that this problem may be fixed 
by equipping DGLAP evolution with additional perturbative BFKL-like resummation of 
small-x effects [32,33], which seem to work in the right direction to reconcile data and 
theory. However, no analysis of data using small- a; resummed DGLAP evolution has been 
performed to date. In any case, it should also be taken into account that DGLAP has 
no predictive power towards small-x, since all the x-dependence of the parton distribution 
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functions is encoded in the initial conditions for the evolution. On the other hand, anal- 
yses relying on the use of fcj-factorization using an unified DGLAP/BFKL resummation 
approach -thus beyond coUinear factorization but still in the linear regime- also provide 
a good description of HERA data for the longitudinal structure function, which is more 
sensitive to saturation effects. Finally, models [34,35] based on the non-pertubative Regge 
calculus have also proven able to describe e+p data on structure functions. 

Thus, and despite the indications for the possible presence of saturation effects outlined 
above, it is not clear from present analyses whether current data signals the breakdown 
of coUinear factorization nor the onset of non-linear corrections to the QCD evolution 
equations. Further, it remains unclear what is the precise kinematic region where these 
two regimes of QCD should be distinguishable. The answer to these questions demands 
either that new data on different observables at smaller values of x is obtained -that is the 
purpose of the proposed experimental facilities as the LHeC [36] or the EIC [37]- or that 
the more accurate data obtained by HI and ZEUS collaborations is analyzed 

Such is precisely the goal of this work. Here we extend our previous analysis by 
including the recent data on reduced cross section in e+p collisions at HERA as given 
by the combined analysis of the HI and ZEUS collaborations [38], which supersedes older 
HERA data on inclusive and longitudinal structure functions. Data from E665 and NMC 
experiments arc kept in the analysis. Besides of extending the explored kinematic range, 
the recent measurements performed at HERA have two clear advantages over previous 
data sets. On one hand, the combination method employed in the analysis reduces the 
systematic uncertainties, resulting in improved accuracy and smaller error bars, thus posing 
more constraining conditions on models. On the other, the fact that data is given directly 
in terms of the reduced cross sections 0"^., which is the experimentally measured quantity, 
eliminates the theoretical bias in the extraction of and F2 from data. 

Besides the analysis of new HERA data, the main novelty in this work is the inclusion 
of dynamical heavy quarks, charm and beauty. It is known from direct experimental 
measurements [39-44] that charm and beauty contribute substantially to the total cross 
section and, therefore, their contribution should not be neglected. Technically, the charm 
and beauty contribution can be obtained by simply extending the sum over quark flavors 
in the dipole model Eq. (2.4). However, as already noticed in [11,45], such straightforward 
inclusion of heavy quarks has strong effects on the fit parameters. In particular it tends 
to strongly reduce the saturation scale of the proton. This is related to the fact that the 
large charm (or beauty) mass acts as an infrared regulator for the fluctuation of the virtual 
photon into a quark-antiquark dipole, thus leaving less room for interaction in the kinematic 
region where unitarity effects are expected to be important. Other fit parameters also 
change noticeably after such direct inclusion of heavy quarks, thus blurring their physical 
interpretation. We notice that a simple modification on the assumptions concerning the 
normalization of the heavy quark contribution or, equivalently, concerning the average 
radius of their transverse distribution, fixes the problem of stability of the fits. For a 
consistent treatment of dynamical heavy quark effects we consider not only its contribution 
to the DIS cross section but we also implement a variable flavor scheme for the beta function 
to properly incorporate such effects in the running of the coupling. 
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This work is structured as follows: Section 2 is devoted to a brief review of our theoret- 
ical setup, which relies on the dipole model formulation of the e+p scattering process and 
in the use of the rcBK equation to describe the small-x dynamics of the dipole scattering 
amplitude. There we discuss the free parameters in the fit, together with our choice of ini- 
tial conditions for the solution of the rcBK equation. The implementation of the variable 
flavor number scheme for the running of the coupling as well as the infrared regularization 
of the coupling are discussed in section 2.2. The experimental data included in the fits and 
the numerical method devised to perform the global fits are discussed in section 3. Our 
results are presented in section 4, where we first present the fits including only light quarks 
in the analysis. We then include the effects of charm and beauty, finding in both cases a 
good description of data. Finally, we wrap up with summary and conclusions. 

2. Setup 

In this section we briefiy review the main ingredients needed for the calculation of the 
inclusive and longitudinal DIS structure functions, which was extensively discussed in our 
previous paper [25]. Neglecting the contribution from Z boson exchange, only relevant 
at much larger than those considered in this work, the reduced cross section can be 
expressed in terms of the inclusive, F2, and longitudinal, Fj,, structure functions: 

X, Q2) = F2(x, Q") - ^^—f-—^Fi^{x, Q^), (2.1) 

where y = Q^/ (s x) is the inelasticity variable and ^/s the center of mass collision energy. 
In turn, at x ^ 1, the inclusive and longitudinal structure functions can be expressed as 

F2{x,Q^)= f {aT + CTL) , (2.2) 

47r^aem 

Fl{x,Q^) = -4 (7L. (2.3) 

Here (Jt,l stands for the virtual photon-proton cross section for transverse (T) and longi- 
tudinal (L) polarization of the virtual photon. In the dipole model, valid at high energies 
or small x, one writes [9, 10]: 

(^t,l{x,Q^) = 2Y, dz J dhdr\^j^{ef,mf,z,Q^r)\''M{h,r,x), (2.4) 

where ^ is the light-cone wave function for a virtual photon to fluctuate into a quark- 
antiquark dipole of quark flavor /. Note that ^j^^ only depends on the quark flavor / 
through the quark mass m/, and electric charge e/ (see e.g. [11] for explicit expressions 
to lowest order in Oem)- A/'(b, r, x) is the imaginary part of the dipole-target scattering 
amplitude, with r the transverse dipole size and b the impact parameter of the collision. 
The study of impact parameter dependence of the dipole amplitude is controlled by long- 
range, non-perturvative phenomena rooted in the physics of confinement and thus is not 
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tractable pcrturbatively. Following other works, we shall assume an average over impact 
parameter through the replacement: 

2 / db 



Jdb^ao, (2.5) 



where (Tq has the meaning of (half) the average transverse area of the quark distribution 
in the transverse plane and will be one of the free parameters in the fit. Attempts to go 
beyond the translational invariant approximation in the BK equation have been recently 
presented in [46,47]. Finally, we shall also include in the fits available data for the charm 
contribution to the inclusive structure function F2c, which can be calculated by considering 
only the charm contribution in Eqs (2.2) and (2.4). Details about the normalization and 
initial conditions for F2c are given below. In order to approach safely the photoproduction 
region, we shall also consider the standard kinematic shift in the definition of Bjorken-a; [11]: 

/ 4m^\ 

x = xh+-^\. (2.6) 

The mass of the three light quarks is taken to be mi = 0.14 GeV in some cases or left as a 
free fit parameter, whereas that of charm and beauty are taken to be rricharm = 1.27 GeV 
and rubeauty = 4.2 GeV respectively [48]. 

2.1 BK equation with running coupling 

The main dynamical input in this work is the rcBK equation, which corresponds to the 
large- ATp limit of the full B-JIMWLK equations. It resums to all orders leading radiative 
corrections in ag ln(l/a;) and also a subset of the full next-to-leading order corrections [21], 
namely running coupling corrections. The impact parameter independent BK equation 
reads 



J drii^™°(r,ri,r2) 



dM{r, x) 
d ln(xo/a;) 

X [J\f{ri,x)+J\f{r2,x)-J\f{r,x)-J\f{ri,x)J\f{r2,x)], (2.7) 
with the evolution kernel including running coupling corrections given by [18] 



K-'^(r,ri,r2) 



27r2 



(2.8) 



where = r — ri and xq is the value of x where the evolution starts. In our case Xq = 0.01 
will be the highest experimental value of x included in the fit. 

2.2 Variable flavor scheme and regularization of the coupling 

The coupling in the rcBK kernel Eq. (2.8) is given, for a given number of active quark 
fiavors n/, by 

a,n,(r^) = —. ^ , (2.9) 
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where 

l3o,nf = n-^nf. (2.10) 

Here, the constant under the logarithm accounts for the uncertainty inherent to the 
Fourier transform from momentum space, where the original calculation of the quark part 
of the /9 function was performed [18,49], to coordinate space. It will be one of the free 
parameters in the fits. 

In both our previous analysis [25] and for the fits in subsection 4.1 only light quarks 
were taken as contributing to the DIS cross section. In this case, only fluctuations of the 
virtual photon wavefunction in Eq. (2.4) into dipoles of light quark flavor were included in 
the calculation. Consistently, only light quark loops should be included in the calculation 
of the running coupling Eq. (2.9). Thus, the number of active flavors in Eq. (2.9) is taken 
to be fixed and equal to the number of light quarks nj = 3. 

Since the rcBK equation is an integro-differential equation where the phase space for 
all dipole sizes is explored, including arbitrarily large dipole sizes (which correspond to 
emission of gluons with arbitrarily small transverse momenta), a prescription to regulate 
the coupling in the infrared is needed. We freeze the coupling to two constant values 
afr = 0.7 and 1 for dipole sizes larger than the scale at which the running coupling reaches 

afr. 

When heavy quark (charm and beauty) contributions are included in the calculation 
of the DIS cross section, as it is the case for the fits in subsection 4.2, fluctuations of the 
virtual photon wavefunction in Eq. (2.4) into dipoles of heavy quark flavor are allowed. 
Accordingly, such contributions should be accounted for in the computation of the running 
coupling Eq. (2.9). Thus, the number of active flavors Uf in Eq. (2.9) should be set to 
the number of quark flavors lighter than the momentum scale associated with the scale 
at which the coupling is evaluated //^ = 4C^/r^. The setup of this variable fiavor 
scheme is completed by matching the branches of the coupling with adjacent n/ at the 
scale corresponding to the quark masses = 4C'^/m^. For the 1-loop accuracy at which 
the coupling Eq. (2.9) is evaluated, the matching condition is simply given by 

as,n/-i(r*) = as,nfirl) , (2.11) 

which results in 

An,-i = {nif)'^^' (A,,)""-/-^ . (2.12) 

The values of the A„^. , A3 in the fixed nj scheme and A3, A4, and A5 for variable nj are 
determined by using an experimentally measured value of as as reference. It is a well known 
fact that the running of the QCD coupling evaluated to 1-loop is of insufficiently accuracy to 
describe the experimental observed coupling evolution. Thus, different choices of reference 
measurement will result in slightly different values for the A„^ . To take into account such 
uncertainty, in some of the fits wc will use as reference point the experimentally measured 
value of as at the mass, whereas in other fits the measured value of the coupling at the 
r mass will be taken as the reference scale. 



-6- 



2.3 Initial conditions for the evolution 



Finally, to complete all the ingredients needed for the calculation of the reduced cross 
section Eq. (2.1) we need to specify the initial conditions for the rcBK evolution equation 
Eq. (2.7). Similarly to our previous work we consider GBW initial conditions, inspired in 
the phenomenological model of [11]: 



AT' 



■GBW 



{r,x = xo) = 1 — exp 



(2.13) 



and MV initial conditions, which originate from a semiclassical calculation of multiple 
rescatterings [50]: 



M''^^{r,x = xo) = 1-exp 



In 



1 
r A 



+ e 



(2.14) 



The physical meaning of the different parameters in Eqs (2.13) and (2.14) is the following: 
Qso is the saturation scale at the largest value of x considered in the analysis, xo = 0.01, 
while 7 is an additional parameter that controls the steepness of the fall-off of the dipole 
amplitude with decreasing r. It should be noted that the factor A under the logarithm 
in the MV initial conditions corresponds to the infrared cutoff of the dipole-nucleon cross 
section at the level of two gluon exchange or in the semiclassical limit. Thus, it does not 
need to be equal to the A„^ in the running of the coupling. However, we opt to set it equal 
to A3. 

In order to further explore the space of initial conditions we shall consider a third family 
of i.e., the scaling i.e. which is generated by the evolution itself. It is a well known result 
that the asymptotic solutions of the rcBK equation are universal, i.e., they are independent 
of the initial conditions [22,51 53]. Moreover, such asymptotic solutions present the feature 
of scaling, i.e. they do no longer depend on two kinematic variables r and Y, but rather on 
a single dimensionless scale, the scaling variable r = r Qs{Y). In other words, the evolution 
generates a universal shape for the dipole amplitude at asymptotically large rapidities 



M{r,Y :^ 1) ^ M'^^'ir = r Qs{Y)) 



(2.15) 



Since the analytic form of the universal shape jV'*'^"' is not known, the implementation of 
the scaling i.e. is done numerically: we solve Eq. (2.7) up to large rapidities, which we 
set to be y = 80. Then the obtained solution is rescaled by the corresponding value of 
the saturation scale, i.e we replace r = rQs{Y) rQso in Eq. (2.15), where Qsq carries 
again the meaning of initial saturation scale at x = xq Thus, the scaling i.e. is essentially 
a one-parameter family of solutions, the only free parameter being the initial saturation 
scale. 



2.4 Parameters for fits with heavy quarks 

As discussed at the beginning of this section, we replace the two-dimensional integral over 
impact parameter in Eq. (2.4) by a dimensionful scale ctq which sets the normalization and 
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can be interpreted as the average transverse size of the proton. However, it is not clear 
a priori whether such average area should be the same for quarks (valence or sea) and 
gluons. Indeed it has been suggested that the glue distribution inside nucleons may be 
located inside hot spots of small radius ~ 0.2 0.3 fm [54]. Also, data on the exponential 
slope of the momentum transfer dependence of exclusive vector meson production (see [55] 
and references therein) provide further support the picture of a smaller effective area for 
gluons than for valence quarks. Here we take as a working hypothesis the possibility that 
the effective transverse size of the heavy quark distribution, which we expect to follow the 
gluon one, may be different to that of light quarks. Accordingly, we introduce two different 
normalization constants for the total cross section, one for charm and beauty, ctq^"''^ and 
other for the three light quarks, ctq: 




As we shall discuss in section 4, such assumption is not only a physically well motivated 
one, but it turns out to be necessary in order to attain a good description of data, and 
also for the stability of the fits with respect to the inclusion or not of the heavy quark 
contribution. Finally, the superscripts light and heavy in the dipole scattering amplitudes 
in Eq. (2.16) refers to the fact that we may consider different initial values of the parameters 
in the initial condition for light and heavy quarks. 

2.5 Summary of the theoretical setup and free parameters 

In summary, we will calculate the reduced cross section and the charm and beauty con- 
tribution to the inclusive structure functions according to the dipole model under the 
translational invariant approximation Eq. (2.16). The small-x dependence is completely 

described by means of the BK equation including running coupling corrections, Eqs. (2.7- 
2.8), for which three different initial conditions GBW, MV and scaling are considered. All 
in all, the free parameters to be fitted to experimental data are: 

• (Jo : The total normalization of the cross section in Eq. (2.16). 

• Q^Q : The saturation scale of the proton at the highest experimental value of Bjorken- 
X included in the fit, xq = lO^^, in Eqs. (2.13) and (2.14). 

• C^: The parameter relating the running of the coupling in momentum space to the 
one in dipole size in Eq. (2.9). 

• 7 : The anomalous dimension of the initial condition for the evolution in Eqs. (2.13) 
and (2.14). 

The fits with heavy quarks introduce additional free parameters, cJq^"^^, Qq*^"^^ and ^^'■^"■'"y ^ 
with physical meaning analogous to that of the corresponding parameters listed above. 
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3. Numerical method and experimental data 

The fit includes data on different observables and from different experiments: 

• Data on the inclusive structure function F2 by the E665 [56] (FNAL), the NMC [57] 
(CERN-SPS) collaborations. 

• Data for the reduced cross section cr^ from the combined analysis of the HI and 
ZEUS collaborations [38]. 

• Data on the charm contribution to the total structure function i*2c [40-44]. Even if 
the beauty contribution to the reduced cross sections is considered, we do not include in 
the fits the few available data on F2b. They have large error bars and influence very little 
the fit output. 

We have considered data for x < 10~^ and for all available values of Q^, 0.045 GeV^ < 
< 50 GeV^. The only published direct measurement of the longitudinal structure func- 
tion Fl{x,Q^) were obtained recently by the HI and ZEUS collaborations [58,59], and it 
is not included in the fit. All in all, 325 data points are included in the fits with only 
light quarks and 329 data points in fits including heavy quarks and the results for F2c- 
It should be noted that data points not fulfilling the condition x < 10~^ see, Eq. (2.6), 
more constraining when charm quark is considered, were excluded from the data set, hence 
the relatively small difference in the total number of data points for fits with or without 
F2c data. Statistical and systematic uncertainties were added in quadrature, and relative 
normalization uncertainties not considered. It was checked in [38] that for the DGLAP 
analysis performed there, this rough treatment of uncertainties did not result in sizable 
differences with respect to the rigorous treatment of statistical, systematic and normaliza- 
tion uncertainties separately which is far more demanding on computer resources. Since 
the minimization algorithms require a large number of calls to the function we have imple- 
mented a parallelization of the numeric code. Finally, the BK evolution equation including 
running coupling corrections is solved using a Runge-Kutta method of second order with 
rapidity step Ahy = 0.05. For further details see [22]. 

4. Results 

4.1 Fits with only light quarks 

We first perform fits to data accounting only for the light quark contribution in Eq. (2.4), 

f = u,d, s. For these fits we take Uf = 3 in the running of the coupling. As the default 
setting we adjust Aqcd to reproduce as{mz) correctly. We assessed the effect of this choice 
by also performing fits where as{mr) was alternatively taken as the reference, obtaining 
similarly good fits. In all these fits we impose an upper cut = 50 in the data. However, 
we note that the fits are stable after extending our analyses to data with larger values 
of < 100 GeV^, albeit at the prize of getting slightly worse x^/d-O-f-- The values of 
the free parameters obtained from the fits to data are presented in Table 1, and a partial 
comparison with the experimental data for the reduced cross section is presented in Fig. 
1 (left plot). Note that, with several combinations resulting in fits of very similar quality. 
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Table 1: Parameters from fits with only light quarks to data with x < 10~^ and for all available 
values of < 50 GeV^ for different initial conditions, fixed values of the coupling in the infrared 
afr = 0.7 and 1 and light quark masses either taken fixed m; = 0.14 GcV or left as a free parameter. 
Fits a', c' and h' correspond to taking the r mass as reference scale for the running of the coupling. 
Units: QIq and mf are in GeV^ and ao in mb. 

both in this subsection and in the following one we will only show in the plots the results 
from some selected fits, not a full survey of them. Several comments are in order: 

First, all the diflFerent fits with MV or GBW initial conditions yield a good /d.o.f < 
1.28, with a best fit /d.o.f = 1.104, labeled h' in table 1, obtained with MV initial 
condition, afr = 1, and as{'mr) as the reference value for the running coupling. The quality 
of the fits is remarkably good provided the tiny error bars in the new data on reduced cross 
sections (error bars are in most cases smaller than the symbols used in the plot). In turn, 
it was not possible to find any good fit to data using the scaling initial condition Eq. (2.15). 
Most likely, this is due to the much faster evolution speed featured by the scaling initial 
conditions, compared to the GBW or MV ones, for which pre-asymptotic effects slow down 
the evolution considerably. Moreover, the MV i.e. tend to systematically yield better fits 
than the GBW one. This can be taken as an indication that the semiclassical resummation 
of multiple scattering underlying the MV formula is indeed a good estimate of the initial 
condition. 

Next, the sensitivity of the fits to non-perturbative aspects of our calculation encoded 
in the parameters afr, C, the reference scale to determine or the light quark masses mi 
(which acts as an effective IR cutoff for the ^'')'*~*'«« wavefunction) is rather small, as shown 
by the little variation of the fit parameters under changes in the latter and on whether 
they arc left as free fit parameters or not. In particular, the value at which the coupling is 
regularized in the infrared, either 1 or 0.7 does not affect much the fit output. Also, when 
the light quark mass is left as a free parameter it tends to acquire a final value very close 
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to the dei one, mi = 140 MeV. This gives us confideE lat the good agreement with 
data is in driven by the small- j; dynamics encoded i ; rcBK equation rather than 
being due i fine tuning of the remaining parameters. ^ ortantly, the fits parameters 
are similar, in ah cases, to those obtained in our previous work [25]. 



Fit with only light quarks 



Fit including heavy quarks 
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Figure 1: Comparison of experimental data for the reduced cross sections (black squares) in 
different bins with our results (red circles). The results in the left plot correspond to a fit with 
only light flavors and GBW initial condition, entry (a) in Table 1. The results in the right plot 
include the contribution of charm and beauty quarks and correspond to fit (a') in Table 2. 



4.2 Inclusion of heavy quarks in the fits. 

In this section we present the fits performed including the contribution of charm and 
beauty quarks into Eq. (2.4). As discussed earlier, we find that in order to obtain a good 
description of data while keeping the stability of the fit parameters for light quarks it 
is necessary to assume that the overall normalization of the heavy quark contribution to 
the reduced cross section is different to the one for light quarks. This translates into the 
introduction of a new free parameter, Uq'^"'^^ , which turns out to be smaller than co, the 
corresponding normalization for light quarks. This can be interpreted as the average radius 
of the heavy quark distribution being smaller than the one for light quarks. In principle, 
there is no reason a priori why such average radius should be the same for charm and 
beauty quarks. On the contrary, one may expect a smaller size of the effective beauty 
distribution on account of its larger mass. This would suggest the introduction of two 
different normalization parameters for charm and beauty ctqc and aob, as well as, maybe, 
different initial conditions for the evolution for each heavy quark flavor. However, the 
paucity of data on or related observables able to independently constrain the free 
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Table 2: Parameters from fits including cliarm and beauty contributions to data with x < 10~ 
and and < 50 GeV^ for different initial conditions and fixed values of the coupling in the infrared 
afr = 0.7 and 0.1. Light quark masses are fixed to mi = 0.14 GeV in some fits and left as a free 
parameter in others. The fit a' corresponds to taking the r mass as reference scale for the running 
of the coupling. The units: (5^o(c) ^^'^ ™f GeV^, while those of (To(c) are mb. 



parameters associated to the beauty quark prevents us of carrying out a more detailed 
characterization of its contribution to the data included in the fit. Thus, we assume that 
the free parameters associated to heavy quarks, including tlic overall normalization, is the 
same for charm and beauty. Wc have checked that such assumption has a very little effect 
on the fit output by completely removing the beauty contribution to F2 and 0"^. However, 
we finally decided to include it in the fits in order to be consistent with the variable flavor 
scheme used for the running of the coupling, which allows the contribution of dynamical b 
quarks to the QCD beta function. 

Our fit results are shown in Table 2, and a comparison with data for ar is shown in 
the right plot Fig 1. We obtain an equally good description of data as with fits with only 
light quarks, as can be seen comparing the left and right plots in Fig 1. However, the 
/d.o.f. < 1.4 are shghtly larger than for the fits with only light quarks. This is maybe 
due to what seems to be a systematic deviation between different data sets on and the 
charm contribution to the reduced cross section arc, as can be observed in Fig 2, where 
we compare our rcsTilts with experimental data. The arguments presented before on the 
stability of the fits with respect to variations in the infrared regulation of the coupling 
or the reference scale to determine Aqqd also hold in the case of fits with heavy quarks. 
On the other hand, when left as a free parameter the mass of the light quark tends to 
acquire a smaller value than it did in the fits with only light quarks. Concerning the 
initial conditions for the evolution, they arc very similar for light and heavy quarks. In 
particular, the corresponding initial saturation scales, Qso and QsOc take on very similar 
values in all fits. However, the steepness of the initial condition encoded in the parameter 
7(c) is systematically larger for light than for heavy quarks for both GBW and MV initial 



- 12 - 



Comparison with data on F^^ and o^, 
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Figure 2: Comparison of experimental data for F2c (black squares) and Crc (red squares) in 
different bins with our results (cyan circles), corresponding to fit (a') in Table 2. 

conditions. 

4.3 Comparison with 

In Fig 3. we present a comparison of our results for the longitudinal structure function 
Fl with the available data at small-x and for different bins. The theoretical results 
were obtained using the dipole parametrizations corresponding to fits (e) and (a) in Tables 
1 and 2 respectively, although we have checked that all the others provide equally good 
comparisons with data. The agreement with data is good, provided the relatively large 
error bars in experimental data. 

5. Conclusions 

In this paper we have presented an analysis of the available data on the several inclusive 
structure functions and reduced cross section measured in e+p collisions at small-x. This 
proves the ability of the rcBK equation, the main dynamical ingredient in our approach, 
to account for the x-dependence of the available data, including the high quality data on 
reduced cross sections provided by the combined analysis of the HI and ZEUS Collabora- 
tions. We thus offer additional indications for the presence of non-linear saturation effects 
in present data and, thereby, sharpen the CGC approach to high-energy QCD scattering as 
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Figure 3: Comparison of experimental data for Fl from the HI (full circles) and ZEUS (triangles) 
collaborations with the theoretical results corresponding to a fit with only light quarks and MV i.e. 
(solid line, labeled (e) in Table 1) and a fit with MV i.e. and including heavy quarks (dashed line, 
labeled (a) in Table 2) . 

a practical phenomenological tool. We have also shown how the inclusion of heavy quarks, 
both at the level of their contribution of the QCD beta function in the running of the 
coupling as well as to the total 7*-proton cross sections can be naturally incorporated in 
the dipole formalism under the assumption of a smaller size of the heavy quark effective 
distribution. The dipole scattering amplitude solving the rcBK equation stemming from 
the parameter sets in Tables 1 and 2 shall be publicly available in the form of numeric 
Fortran routines at the website http://www-fp.usc.es/phenom/software.html. 
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